# =============================================================================
# May 6 2025
# R code used to produce all results, tables, and figures in the manuscript
# Rebecca Cordell
# Unpacking the Role of In-Group Bias in US Public Opinion on Human Rights Violations
# American Journal of Political Science
# https://doi.org/10.7910/DVN/TGAL7M
# =============================================================================

# Clear work environment
rm(list=ls())

# Install Packages
#install.packages("dplyr")
#install.packages("ggplot2")
#install.packages("forcats")
#install.packages("ggpubr")
#install.packages("ggeasy")
#install.packages("lemon")
#install.packages("sandwich")
#install.packages("survey")
#install.packages("lmtest")
#install.packages("remotes")
#remotes::install_version("cregg", version = "0.4.0")

# Required Packages
library("dplyr")
library("ggplot2")
library("forcats")
library("ggpubr")
library("sandwich")
library("survey")
library("lmtest")
library("cregg")
library("ggeasy")
library("lemon")

options(scipen = 999)
options(warn=-1)

# -----------------------------------------------------------------------------
# Read in data
# -----------------------------------------------------------------------------

hr_survey<-read.csv("cordell_ingroupbiashumanrights_data.csv", header=TRUE, stringsAsFactors = FALSE)

# =============================================================================
# Figure 2: Effect of group identity variables on respondents’ disapproval
# =============================================================================

# Convert regression variables into factors
factor_vars <- c("perp_match", "agent", "type", "scope", "targ_nonstate",
                 "targ_race_match", "targ_relig_match", "targ_citiz_match",
                 "frame", "elite_match")
hr_survey[factor_vars] <- lapply(hr_survey[factor_vars], as.factor)

# -----------------------------------------------------------------------------
# Calculate marginal means
# -----------------------------------------------------------------------------

# Model 1
fig2_m1 <- cregg::cj(hr_survey, outcome1 ~ perp_match + agent + type + scope + targ_nonstate + targ_race_match + targ_relig_match + targ_citiz_match + frame + elite_match, id = ~id, estimate = "mm")

# Subset group identity dummy variables
group_vars <- c("perp_match", "targ_race_match", 
                "targ_relig_match", "targ_citiz_match", 
                "elite_match")
fig2_m1 <- fig2_m1[fig2_m1$feature %in% group_vars,]

# Model 2
fig2_m2 <- cregg::cj(hr_survey, outcome2 ~ perp_match + agent + type + scope + targ_nonstate + targ_race_match + targ_relig_match + targ_citiz_match + frame + elite_match, id = ~id, estimate = "mm")

# Subset group identity dummy variables
fig2_m2 <- fig2_m2[fig2_m2$feature %in% group_vars,]

# -----------------------------------------------------------------------------
# Prepare figures
# -----------------------------------------------------------------------------

# Combine models
fig2_all <- rbind(fig2_m1, fig2_m2)

# Create group identity variable labels
variable_labels <- c(
  "perp_match" = "Perpetrator (partisanship)",
  "targ_race_match" = "Target (race)",
  "targ_relig_match" = "Target (religion)",
  "targ_citiz_match" = "Target (citizenship)",
  "elite_match" = "Elite cue (partisanship)"
)
fig2_all <- dplyr::bind_rows(
  fig2_all,
  do.call(rbind, lapply(names(variable_labels), function(variable) {
    data.frame(
      outcome = c("outcome1", "outcome2"),
      variable = variable_labels[variable],
      level = variable_labels[variable],
      estimate = NA, lower = NA, upper = NA
    )
  }))
)

# Create respondents variable (in-groups vs. out-groups)
fig2_all <- fig2_all %>%
  mutate(
    Respondents = case_when(
      grepl("_in", level) ~ "In-group",
      grepl("_out", level) ~ "Out-group",
      level %in% c("Perpetrator (partisanship)", "Target (race)", 
                   "Target (religion)", "Target (citizenship)", "Elite cue (partisanship)") ~ "In-group",
      TRUE ~ level
    ),
    Respondents = factor(Respondents, levels = c("In-group", "Out-group"))
  )

# Set factor order
level_order <- c(
  "Perpetrator (partisanship)", "perp_in", "perp_out",
  "Target (race)", "race_in", "race_out",
  "Target (religion)", "relig_in", "relig_out",
  "Target (citizenship)", "citiz_in", "citiz_out",
  "Elite cue (partisanship)", "elite_in", "elite_out"
)
fig2_all <- fig2_all %>%
  mutate(level = factor(level, levels = level_order) %>% fct_rev())

# Separate models
fig2_m1 <- filter(fig2_all, outcome == "outcome1")
fig2_m2 <- filter(fig2_all, outcome == "outcome2")

# -----------------------------------------------------------------------------
# Create figures
# -----------------------------------------------------------------------------

# Create panel (a) Disapproval forced-choice
fig2_m1_plot <- ggplot(fig2_m1, aes(x = level, y = estimate, colour = Respondents)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), size = 0.3, width = 0.2) +
  scale_x_discrete(labels = c(
    "", "Elite cue (partisanship)", "", 
    "", "Target (citizenship)", "", 
    "", "Target (religion)", "", 
    "", "Target (race)", "", 
    "","Perpetrator (partisanship)",""
  )) +
  xlab('') + ylab('Marginal mean') +
  coord_flip() +
  geom_hline(yintercept = 0.5, size = 0.2, color = "black") +
  theme_classic(base_size = 10) +
  scale_colour_grey() +
  easy_center_title() +
  ggtitle("(a) Disapproval forced-choice") +
  theme(
    axis.text.x = element_text(colour = "black"),
    axis.text.y = element_text(colour = "black"),
    axis.title.x = element_text(margin = margin(t = 10)),
    axis.ticks.y = element_blank()
  ) +
  labs(color = NULL)

# Create panel (b) Disapproval ratings-based
fig2_m2_plot <- ggplot(fig2_m2, aes(x = level, y = estimate, colour = Respondents)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), size = 0.3, width = 0.2) +
  scale_x_discrete(labels = c(
    "", "Elite cue (partisanship)", "", 
    "", "Target (citizenship)", "", 
    "", "Target (religion)", "", 
    "", "Target (race)", "", 
    "","Perpetrator (partisanship)",""
  )) +
  xlab('') + ylab('Marginal mean') +
  coord_flip() +
  theme_classic(base_size = 10) +
  scale_colour_grey() +
  easy_center_title() +
  ggtitle("(b) Disapproval ratings-based") +
  theme(
    axis.text.x = element_text(colour = "black"),
    axis.text.y = element_text(colour = "black"),
    axis.title.x = element_text(margin = margin(t = 10)),
    axis.ticks.y = element_blank()
  ) +
  labs(color = NULL)

# Print Figure 2
pdf('figure2.pdf', width = 8.26, height = 5.82)
ggarrange(fig2_m1_plot, NULL, fig2_m2_plot,
          nrow = 1, common.legend = T, legend = "bottom", ncol=3, widths = c(1, 0.09, 1))
dev.off()

# -----------------------------------------------------------------------------
# In-text Results
# -----------------------------------------------------------------------------

# "For the forced-choice outcome, abuses by out-group perpetrators had a 6 percentage point higher marginal mean of disapproval (0.53) than those by in-group perpetrators (0.47), a difference three times larger than any other group identity variable."
# 6 percentage point higher
round((fig2_m1$estimate[2]-fig2_m1$estimate[1])*100, digit=0) 
# (0.53)
round(fig2_m1$estimate[2], digit=2)
# (0.47)
round(fig2_m1$estimate[1], digit=2)
# three times larger
race<-round((fig2_m1$estimate[3]-fig2_m1$estimate[4]), digit=2)
relig<-round((fig2_m1$estimate[5]-fig2_m1$estimate[6]), digit=2)
citiz<-round((fig2_m1$estimate[7]-fig2_m1$estimate[8]), digit=2)
elite<-round((fig2_m1$estimate[9]-fig2_m1$estimate[10]), digit=2)
perp<-round((fig2_m1$estimate[2]-fig2_m1$estimate[1]), digit=2)
perp/relig

# "This pattern is consistent in the ratings-based outcome, where out-group perpetrators had a higher marginal mean of disapproval (0.72) than in-group perpetrators (0.70)."
# (0.72) 
round(fig2_m2$estimate[2], digit=2)
# (0.70)
round(fig2_m2$estimate[1], digit=2)

# "Religion produced a 2 percentage point higher marginal mean of disapproval (0.51) for religious in-groups in the forced-choice outcome compared to out-groups (0.49) and a 1 percentage point higher marginal mean of disapproval in the ratings-based outcome (0.69 vs. 0.68)."
# 2 percentage point higher
round((fig2_m1$estimate[5]-fig2_m1$estimate[6])*100, digit=0) 
# (0.51)
round(fig2_m1$estimate[5], digit=2)
# (0.49)
round(fig2_m1$estimate[6], digit=2)
# 1 percentage point higher
round((fig2_m2$estimate[5]-fig2_m2$estimate[6])*100, digit=0)
# (0.69 vs. 0.68)
round(fig2_m2$estimate[5], digit=2)
round(fig2_m2$estimate[6], digit=2)

# =============================================================================
# Figure 3: Effect of Perpetrator (Partisanship) Attribute on Respondents' Disapproval Conditional on Respondents’ Party Affiliation
# =============================================================================

# Convert regression variables into factors
factor_vars <- c("perp", "agent", "type", "scope", "targ_nonstate", "targ_race",
                 "targ_relig", "targ_citiz", "frame", "elite", "resp_dem", "resp_rep")
hr_survey[factor_vars] <- lapply(hr_survey[factor_vars], as.factor)

# -----------------------------------------------------------------------------
# Calculate subgroup marginal means
# -----------------------------------------------------------------------------

# Model 1 
fig3_m1_dem <- cregg::cj(hr_survey, outcome1 ~ perp + agent + type + scope + targ_nonstate + targ_race + targ_relig + targ_citiz + frame + elite, id = ~id, estimate = "mm", by = ~resp_dem)
fig3_m1_rep <- cregg::cj(hr_survey, outcome1 ~ perp + agent + type + scope + targ_nonstate + targ_race + targ_relig + targ_citiz + frame + elite, id = ~id, estimate = "mm", by = ~resp_rep)

# Model 2 
fig3_m2_dem <- cregg::cj(hr_survey, outcome2 ~ perp + agent + type + scope + targ_nonstate + targ_race + targ_relig + targ_citiz + frame + elite, id = ~id, estimate = "mm", by = ~resp_dem)
fig3_m2_rep <- cregg::cj(hr_survey, outcome2 ~ perp + agent + type + scope + targ_nonstate + targ_race + targ_relig + targ_citiz + frame + elite, id = ~id, estimate = "mm", by = ~resp_rep)

# -----------------------------------------------------------------------------
# Prepare Figures
# -----------------------------------------------------------------------------

# Subset perpetrator attribute
filter_principal <- function(df, party_col) {
  df %>%
    filter(feature == "perp", !!sym(party_col) == 1) %>%
    mutate(attrib = "Perpetrator (Partisanship)")
}
fig3_m1_dem <- filter_principal(fig3_m1_dem, "resp_dem")
fig3_m2_dem <- filter_principal(fig3_m2_dem, "resp_dem")
fig3_m1_rep <- filter_principal(fig3_m1_rep, "resp_rep")
fig3_m2_rep <- filter_principal(fig3_m2_rep, "resp_rep")

# Combine models
fig3_dem_all <- rbind(fig3_m1_dem, fig3_m2_dem)
fig3_rep_all <- rbind(fig3_m1_rep, fig3_m2_rep)

# Create respondents variable (in-groups vs. out-groups)
fig3_dem_all <- fig3_dem_all %>%
  mutate(Respondents = case_when(
    grepl("Democrat", level) ~ "In-group",
    grepl("Republican", level) ~ "Out-group",
    level == "Perpetrator (Partisanship)" ~ "In-group",
    TRUE ~ level
  )) %>%
  mutate(Respondents = factor(Respondents, levels = c("In-group", "Out-group")))
fig3_rep_all <- fig3_rep_all %>%
  mutate(Respondents = case_when(
    grepl("Republican", level) ~ "In-group",
    grepl("Democrat", level) ~ "Out-group",
    level == "Perpetrator (Partisanship)" ~ "In-group",
    TRUE ~ level
  )) %>%
  mutate(Respondents = factor(Respondents, levels = c("In-group", "Out-group")))

# Set factor order
fig3_dem_all <- fig3_dem_all %>%
  mutate(level = factor(level, levels = c("Perpetrator (Partisanship)", "A Democrat governor", "The Democrat president", "A Republican governor", "The Republican president")) %>% fct_rev())
fig3_rep_all <- fig3_rep_all %>%
  mutate(level = factor(level, levels = c("Perpetrator (Partisanship)", "A Republican governor", "The Republican president", "A Democrat governor", "The Democrat president")) %>% fct_rev())

# Separate models
fig3_m1_dem <- filter(fig3_dem_all, outcome == "outcome1")
fig3_m2_dem <- filter(fig3_dem_all, outcome == "outcome2")
fig3_m1_rep <- filter(fig3_rep_all, outcome == "outcome1")
fig3_m2_rep <- filter(fig3_rep_all, outcome == "outcome2")

# -----------------------------------------------------------------------------
# Create Figures
# -----------------------------------------------------------------------------

# Create panel (a) Democrat respondents - disapproval forced-choice
fig3_m1_dem_plot <- ggplot(fig3_m1_dem, aes(x = level, y = estimate, colour = Respondents)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), size = 0.3, width = 0.2) +
  scale_x_discrete(labels = c("Republican president", "Republican governor", "Democrat president", "Democrat governor")) +
  xlab('') + ylab('Marginal mean') +
  coord_flip(ylim = c(0.42, 0.58)) +
  geom_hline(yintercept = 0.5, size = 0.2, color = "black") +
  theme_classic(base_size = 10) +
  scale_colour_grey() +
  easy_center_title() +
  ggtitle("(a) Disapproval forced-choice") +
  theme(
    axis.text.x = element_text(colour = "black"),
    axis.text.y = element_text(colour = "black"),
    axis.title.x = element_text(margin = margin(t = 10)),
    axis.ticks.y = element_blank(),
    legend.position = "none"
  ) +
  labs(color = NULL)

# Create panel (b) Democrat respondents - disapproval ratings-based
fig3_m2_dem_plot <- ggplot(fig3_m2_dem, aes(x = level, y = estimate, colour = Respondents)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), size = 0.3, width = 0.2) +
  scale_x_discrete(labels = c("Republican president", "Republican governor", "Democrat president", "Democrat governor")) +
  xlab('') + ylab('Marginal mean') +
  coord_flip(ylim = c(0.6, 0.78)) +
  theme_classic(base_size = 10) +
  scale_colour_grey() +
  easy_center_title() +
  ggtitle("(b) Disapproval ratings-based") +
  theme(
    axis.text.x = element_text(colour = "black"),
    axis.text.y = element_text(colour = "black"),
    axis.title.x = element_text(margin = margin(t = 10)),
    axis.ticks.y = element_blank(),
    legend.position = "none"
  ) +
  labs(color = NULL)

# Create panel (c) Republican respondents - disapproval forced choice
fig3_m1_rep_plot <- ggplot(fig3_m1_rep, aes(x = level, y = estimate, colour = Respondents)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), size = 0.3, width = 0.2) +
  scale_x_discrete(labels = c("Democrat president", "Democrat governor", "Republican president", "Republican governor")) +
  xlab('') + ylab('Marginal mean') +
  coord_flip(ylim = c(0.42, 0.58)) +
  geom_hline(yintercept = 0.5, size = 0.2, color = "black") +
  theme_classic(base_size = 10) +
  scale_colour_grey() +
  easy_center_title() +
  ggtitle("(c) Disapproval forced-choice") +
  theme(
    axis.text.x = element_text(colour = "black"),
    axis.text.y = element_text(colour = "black"),
    axis.title.x = element_text(margin = margin(t = 10)),
    axis.ticks.y = element_blank()
  ) +
  labs(color = NULL)

# Create panel (d) Republican Respondents - disapproval ratings-based
fig3_m2_rep_plot <- ggplot(fig3_m2_rep, aes(x = level, y = estimate, colour = Respondents)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), size = 0.3, width = 0.2) +
  scale_x_discrete(labels = c("Democrat president", "Democrat governor", "Republican president", "Republican governor")) +
  xlab('') + ylab('Marginal mean') +
  coord_flip(ylim = c(0.6, 0.78)) +
  theme_classic(base_size = 10) +
  scale_colour_grey() +
  easy_center_title() +
  ggtitle("(d) Disapproval ratings-based") +
  theme(
    axis.text.x = element_text(colour = "black"),
    axis.text.y = element_text(colour = "black"),
    axis.title.x = element_text(margin = margin(t = 10)),
    axis.ticks.y = element_blank()
  ) +
  labs(color = NULL)

# Print Figure 3
fig3_dem_plot<-ggarrange(fig3_m1_dem_plot, NULL, fig3_m2_dem_plot,
                         nrow = 1, common.legend = F, legend = NULL, ncol = 3, widths = c(1, 0.09, 1))
fig3_dem_plot<-annotate_figure(fig3_dem_plot, top = text_grob("Democrats", size = 14))

fig3_rep_plot<-ggarrange(fig3_m1_rep_plot, NULL, fig3_m2_rep_plot,
                           nrow = 1, common.legend = T, legend = "bottom", ncol = 3, widths = c(1, 0.09, 1))
fig3_rep_plot<-annotate_figure(fig3_rep_plot, top = text_grob("Republicans", size = 14))
pdf('figure3.pdf', width = 8.26, height = 5.82)
ggarrange(fig3_dem_plot, fig3_rep_plot,
          nrow = 2, common.legend = T, legend = "bottom", ncol=1, heights = c(1, 1.1))
dev.off()

# -----------------------------------------------------------------------------
# In-text Results
# -----------------------------------------------------------------------------

# "Among Democrats, forced-choice disapproval (panel a) increased by 4–6 percentage points when violations were perpetrated by an out-group Republican president or governor (0.52, 0.53 respectively) compared to an in-group Democrat president or governor (0.47, 0.48 respectively)."
# 4–6 percentage points
dem_gov_rep_gov<-round((fig3_m1_dem$estimate[2]-fig3_m1_dem$estimate[1])*100, digit=0)
dem_gov_rep_pres<-round((fig3_m1_dem$estimate[4]-fig3_m1_dem$estimate[1])*100, digit=0) 
dem_pres_rep_gov<-round((fig3_m1_dem$estimate[2]-fig3_m1_dem$estimate[3])*100, digit=0) 
dem_pres_rep_pres<-round((fig3_m1_dem$estimate[4]-fig3_m1_dem$estimate[3])*100, digit=0) 
min(dem_gov_rep_gov, dem_gov_rep_pres, dem_pres_rep_gov, dem_pres_rep_pres)
max(dem_gov_rep_gov, dem_gov_rep_pres, dem_pres_rep_gov, dem_pres_rep_pres)
# (0.52, 0.53 respectively)
round(fig3_m1_dem$estimate[4], digit=2)
round(fig3_m1_dem$estimate[2], digit=2)
# (0.47, 0.48 respectively)
round(fig3_m1_dem$estimate[3], digit=2)
round(fig3_m1_dem$estimate[1], digit=2)

# "Among Republicans, the difference is greater: disapproval (panel c) increased by 7–10 percentage points when the perpetrator was an out-group Democrat president or governor (0.54, 0.55 respectively) compared to an in-group Republican president or governor (0.44, 0.47 respectively)."
# 7–10 percentage points
rep_gov_dem_gov<-round((fig3_m1_rep$estimate[1]-fig3_m1_rep$estimate[2])*100, digit=0) 
rep_gov_dem_pres<-round((fig3_m1_rep$estimate[3]-fig3_m1_rep$estimate[2])*100, digit=0) 
rep_pres_dem_gov<-round((fig3_m1_rep$estimate[1]-fig3_m1_rep$estimate[4])*100, digit=0) 
rep_pres_dem_pres<-round((fig3_m1_rep$estimate[3]-fig3_m1_rep$estimate[4])*100, digit=0) 
min(rep_gov_dem_gov, rep_gov_dem_pres, rep_pres_dem_gov, rep_pres_dem_pres)
max(rep_gov_dem_gov, rep_gov_dem_pres, rep_pres_dem_gov, rep_pres_dem_pres)
# (0.54, 0.55 respectively)
round(fig3_m1_rep$estimate[3], digit=2)
round(fig3_m1_rep$estimate[1], digit=2)
# (0.44, 0.47 respectively)
round(fig3_m1_rep$estimate[4], digit=2)
round(fig3_m1_rep$estimate[2], digit=2)

# =============================================================================
# Figure 4: Interaction effect of perpetrator (partisanship) with target and elite cue giver identities on respondents’ disapproval
# =============================================================================

# Convert regression variables into factors
factor_vars <- c("agent", "type", "scope", "targ_nonstate",
                 "targ_race_match", "targ_relig_match", "targ_citiz_match",
                 "frame", "elite_match", "perp_targ_race", "perp_targ_relig", 
                 "perp_targ_citiz", "perp_elite","perp_match")
hr_survey[factor_vars] <- lapply(hr_survey[factor_vars], as.factor)

# -----------------------------------------------------------------------------
# Calculate marginal means interactions
# -----------------------------------------------------------------------------

# Model 1
fig4_m1 <- cregg::cj(hr_survey, outcome1 ~ agent + type + scope + targ_nonstate + targ_race_match + targ_relig_match + targ_citiz_match + frame + elite_match + perp_targ_race + perp_targ_relig + perp_targ_citiz + perp_elite, id = ~id, estimate = "mm", by = ~perp_match)

# Model 2
fig4_m2 <- cregg::cj(hr_survey, outcome2 ~ agent + type + scope + targ_nonstate + targ_race_match + targ_relig_match + targ_citiz_match + frame + elite_match + perp_targ_race + perp_targ_relig + perp_targ_citiz + perp_elite, id = ~id, estimate = "mm", by = ~perp_match)

# Subset group identity dummy variables
interaction_terms <- c("perp_targ_race", "perp_targ_relig", "perp_targ_citiz", "perp_elite")
fig4_m1 <- fig4_m1[fig4_m1$feature %in% interaction_terms,]
fig4_m2 <- fig4_m2[fig4_m2$feature %in% interaction_terms,]

# Exclude NA values
fig4_m1 <- fig4_m1[!is.na(fig4_m1$estimate),]
fig4_m2 <- fig4_m2[!is.na(fig4_m2$estimate),]

# -----------------------------------------------------------------------------
# Prepare Figures
# -----------------------------------------------------------------------------

# Combine models
fig4_all <- rbind(fig4_m1, fig4_m2) 

# Create group identity variable labels
variable_labels <- c("targ_race_match" = "Target (race)",
                     "targ_relig_match" = "Target (religion)",
                     "targ_citiz_match" = "Target (citizenship)",
                     "elite_match" = "Elite cue (partisanship)"
)
fig4_all <- dplyr::bind_rows(
  fig4_all,
  do.call(rbind, lapply(names(variable_labels), function(variable) {
    data.frame(
      outcome = c("outcome1", "outcome2"),
      variable = variable_labels[variable],
      level = variable_labels[variable],
      estimate = NA, lower = NA, upper = NA
    )
  }))
)

# Set factor order
fig4_all_perp_in <- filter(fig4_all, !grepl("perp_out", level))
fig4_all_perp_in <- fig4_all_perp_in %>%
  mutate(level = factor(level, levels = c("Target (race)", "perp_in_race_in", "perp_in_race_out",
                                          "Target (religion)", "perp_in_relig_in", "perp_in_relig_out",
                                          "Target (citizenship)", "perp_in_citiz_in", "perp_in_citiz_out",
                                          "Elite cue (partisanship)", "perp_in_elite_in", "perp_in_elite_out")) %>% fct_rev())
fig4_all_perp_out <- filter(fig4_all, !grepl("perp_in", level))
fig4_all_perp_out <- fig4_all_perp_out %>%
  mutate(level = factor(level, levels = c("Target (race)", "perp_out_race_in", "perp_out_race_out",
                                          "Target (religion)", "perp_out_relig_in", "perp_out_relig_out",
                                          "Target (citizenship)", "perp_out_citiz_in", "perp_out_citiz_out",
                                          "Elite cue (partisanship)", "perp_out_elite_in", "perp_out_elite_out")) %>% fct_rev())

# Create respondents variable (in-groups vs. out-groups)
fig4_all_perp_in <- fig4_all_perp_in %>%
  mutate(Respondents = case_when(
    grepl("_in$", level) ~ "In-group",   
    grepl("_out$", level) ~ "Out-group", 
    level %in% c("Target (race)", "Target (religion)", 
                 "Target (citizenship)", "Elite cue (partisanship)") ~ "In-group", 
    TRUE ~ level
  )) %>%
  mutate(Respondents = factor(Respondents, levels = c("In-group", "Out-group")))
fig4_all_perp_out <- fig4_all_perp_out %>%
  mutate(Respondents = case_when(
    grepl("_out$", level) ~ "Out-group",  
    grepl("_in$", level) ~ "In-group",  
    level %in% c("Target (race)", "Target (religion)", 
                 "Target (citizenship)", "Elite cue (partisanship)") ~ "In-group", 
    TRUE ~ level
  )) %>%
  mutate(Respondents = factor(Respondents, levels = c("In-group", "Out-group")))

# Separate models
fig4_m1_perp_in <- filter(fig4_all_perp_in, outcome == "outcome1")
fig4_m2_perp_in <- filter(fig4_all_perp_in, outcome == "outcome2")
fig4_m1_perp_out <- filter(fig4_all_perp_out, outcome == "outcome1")
fig4_m2_perp_out <- filter(fig4_all_perp_out, outcome == "outcome2")

# -----------------------------------------------------------------------------
# Create Figures
# -----------------------------------------------------------------------------

# Create panel (a) Perpetrator in-group - disapproval forced-choice
fig4_m1_perp_in_plot <- ggplot(fig4_m1_perp_in, aes(x = level, y = estimate, colour = Respondents)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), size = 0.3, width = 0.2) +
  scale_x_discrete(labels = c(
    "", "Elite cue (partisanship)", "", 
    "", "Target (citizenship)", "", 
    "", "Target (religion)", "", 
    "", "Target (race)", ""
  )) +
  xlab('') + ylab('Marginal mean') +
  coord_flip() +
  geom_hline(yintercept = 0.5, size = 0.2, color = "black") +
  scale_y_symmetric(mid = 0.5) +
  theme_classic(base_size = 10) +
  scale_colour_grey() +
  easy_center_title() +
  ggtitle("(a) Disapproval forced-choice") +
  theme(
    axis.text.x = element_text(colour = "black"),
    axis.text.y = element_text(colour = "black"),
    axis.title.x = element_text(margin = margin(t = 10)),
    axis.ticks.y = element_blank(),
    legend.position = "none"
  ) +
  labs(color = NULL)

# Create panel (b) Perpetrator in-group - disapproval ratings-based
fig4_m2_perp_in_plot <- ggplot(fig4_m2_perp_in, aes(x = level, y = estimate, colour = Respondents)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), size = 0.3, width = 0.2) +
  scale_x_discrete(labels = c(
    "", "Elite cue (partisanship)", "", 
    "", "Target (citizenship)", "", 
    "", "Target (religion)", "", 
    "", "Target (race)", ""
  )) +
  xlab('') + ylab('Marginal mean') +
  coord_flip(ylim = c(0.65, 0.75)) +
  theme_classic(base_size = 10) +
  scale_colour_grey() +
  easy_center_title() +
  ggtitle("(b) Disapproval ratings-based") +
  theme(
    axis.text.x = element_text(colour = "black"),
    axis.text.y = element_text(colour = "black"),
    axis.title.x = element_text(margin = margin(t = 10)),
    axis.ticks.y = element_blank(),
    legend.position = "none"
  ) +
  labs(color = NULL)

# Create panel (c) Perpetrator out-group - disapproval forced choice
fig4_m1_perp_out_plot <- ggplot(fig4_m1_perp_out, aes(x = level, y = estimate, colour = Respondents)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), size = 0.3, width = 0.2) +
  scale_x_discrete(labels = c(
    "", "Elite cue (partisanship)", "", 
    "", "Target (citizenship)", "", 
    "", "Target (religion)", "", 
    "", "Target (race)", ""
  )) +
  xlab('') + ylab('Marginal mean') +
  coord_flip() +
  geom_hline(yintercept = 0.5, size = 0.2, color = "black") +
  scale_y_symmetric(mid = 0.5) +
  theme_classic(base_size = 10) +
  scale_colour_grey() +
  easy_center_title() +
  ggtitle("(c) Disapproval forced-choice") +
  theme(
    axis.text.x = element_text(colour = "black"),
    axis.text.y = element_text(colour = "black"),
    axis.title.x = element_text(margin = margin(t = 10)),
    axis.ticks.y = element_blank()
  ) +
  labs(color = NULL)

# Create panel (d) Perpetrator out-group - disapproval ratings-based
fig4_m2_perp_out_plot <- ggplot(fig4_m2_perp_out, aes(x = level, y = estimate, colour = Respondents)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), size = 0.3, width = 0.2) +
  scale_x_discrete(labels = c(
    "", "Elite cue (partisanship)", "", 
    "", "Target (citizenship)", "", 
    "", "Target (religion)", "", 
    "", "Target (race)", ""
  )) +
  xlab('') + ylab('Marginal mean') +
  coord_flip(ylim = c(0.65, 0.75)) +
  theme_classic(base_size = 10) +
  scale_colour_grey() +
  easy_center_title() +
  ggtitle("(d) Disapproval ratings-based") +
  theme(
    axis.text.x = element_text(colour = "black"),
    axis.text.y = element_text(colour = "black"),
    axis.title.x = element_text(margin = margin(t = 10)),
    axis.ticks.y = element_blank()
  ) +
  labs(color = NULL)

# Print Figure 4
fig4_perp_in_plot<-ggarrange(fig4_m1_perp_in_plot, NULL, fig4_m2_perp_in_plot,
                             nrow = 1, common.legend = F, legend = NULL, ncol = 3, widths = c(1, 0.09, 1))
fig4_perp_in_plot<-annotate_figure(fig4_perp_in_plot, top = text_grob("Perpetrator in-group", size = 14))

fig4_perp_out_plot<-ggarrange(fig4_m1_perp_out_plot, NULL, fig4_m2_perp_out_plot,
                              nrow = 1, common.legend = T, legend = "bottom", ncol = 3, widths = c(1, 0.09, 1))
fig4_perp_out_plot<-annotate_figure(fig4_perp_out_plot, top = text_grob("Perpetrator out-group", size = 14))
pdf('figure4.pdf', width = 8.26, height = 5.82)
ggarrange(fig4_perp_in_plot, fig4_perp_out_plot,
          nrow = 2, common.legend = T, legend = "bottom", ncol=1, heights = c(1, 1.1))
dev.off()

# -----------------------------------------------------------------------------
# In-text Results
# -----------------------------------------------------------------------------

# "For violations targeting a respondents’ racial in-group, forced-choice disapproval was 6 percentage points lower for in-group perpetrators (0.49) compared to out-group partisan perpetrators (0.55)."
# 6 percentage points lower
round((fig4_m1_perp_in$estimate[1]-fig4_m1_perp_out$estimate[1])*100, digit=0) 
# (0.49)
round(fig4_m1_perp_in$estimate[1], digit=2)
# (0.55)
round(fig4_m1_perp_out$estimate[1], digit=2)

# "However, in-group perpetrators targeting religious or citizenship in-groups faced little backlash, with 7–8 percentage points lower forced-choice disapproval (both 0.47) than out-group perpetrators (0.54, 0.55 respectively)."
# 7–8% lower
round((fig4_m1_perp_in$estimate[5]-fig4_m1_perp_out$estimate[5])*100, digit=0) 
round((fig4_m1_perp_in$estimate[3]-fig4_m1_perp_out$estimate[3])*100, digit=0) 
# (both 0.47)
round(fig4_m1_perp_in$estimate[3], digit=2)
round(fig4_m1_perp_in$estimate[5], digit=2)
# (0.54, 0.55 respectively)
round(fig4_m1_perp_out$estimate[5], digit=2)
round(fig4_m1_perp_out$estimate[3], digit=2)

# =============================================================================
# Figure 5: Additive effect of multiple target and elite cue giver identities on respondents’ disapproval
# =============================================================================

# Convert regression variables into factors
factor_vars <- c("perp_match", "agent", "type", "scope", "targ_nonstate", "targ_race_match", 
                 "targ_relig_match", "targ_citiz_match", "frame", "elite_match", 
                 "shared_ident")
hr_survey[factor_vars] <- lapply(hr_survey[factor_vars], as.factor)

# Create shared identities categorical variable
hr_survey<-hr_survey[hr_survey$shared_ident!=4,]
hr_survey<-hr_survey[!is.na(hr_survey$shared_ident),]
hr_survey$shared_ident <- factor(hr_survey$shared_ident, 
                                      levels = 0:3, 
                                      labels = c("shared_0", "shared_1", "shared_2", "shared_3"))

# -----------------------------------------------------------------------------
# Calculate marginal means
# -----------------------------------------------------------------------------

# Model 1
fig5_m1 <- cregg::cj(hr_survey, outcome1 ~ perp_match + agent + type + scope + targ_nonstate + targ_race_match + targ_relig_match + targ_citiz_match + frame + elite_match + shared_ident, id = ~id, estimate = "mm")

# Model 2
fig5_m2 <- cregg::cj(hr_survey, outcome2 ~ perp_match + agent + type + scope + targ_nonstate + targ_race_match + targ_relig_match + targ_citiz_match + frame + elite_match + shared_ident, id = ~id, estimate = "mm")

# Subset shared identity variable
group_vars <- c("shared_ident")
fig5_m1 <- fig5_m1[fig5_m1$feature %in% group_vars,]
fig5_m2 <- fig5_m2[fig5_m2$feature %in% group_vars,]

# -----------------------------------------------------------------------------
# Prepare figures
# -----------------------------------------------------------------------------

# Combine models
fig5_all <- rbind(fig5_m1, fig5_m2)

# Set factor order
level_order <- c(
  "shared_3", "shared_2", "shared_1", "shared_0"
)
fig5_all <- fig5_all %>%
  mutate(level = factor(level, levels = level_order) %>% fct_rev())

# Separate models
fig5_m1 <- filter(fig5_all, outcome == "outcome1")
fig5_m2 <- filter(fig5_all, outcome == "outcome2")

# -----------------------------------------------------------------------------
# Create figures
# -----------------------------------------------------------------------------

# Create panel (a) Disapproval forced-choice
fig5_m1_plot <- ggplot(fig5_m1, aes(x = level, y = estimate)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), size = 0.3, width = 0.2) +
  scale_x_discrete(labels = c("None", "One", "Two", "Three")) +
  xlab('') + ylab('Marginal mean') +
  coord_flip() +
  geom_hline(yintercept = 0.5, size = 0.2, color = "black") +
  scale_y_symmetric(mid = 0.5) +
  theme_classic(base_size = 10) +
  scale_colour_grey() +
  easy_center_title() +
  ggtitle("(a) Disapproval forced-choice") +
  theme(
    axis.text.x = element_text(colour = "black"),
    axis.text.y = element_text(colour = "black"),
    axis.title.x = element_text(margin = margin(t = 10)),
    axis.ticks.y = element_blank()
  )

# Create panel (b) Disapproval ratings-based
fig5_m2_plot <- ggplot(fig5_m2, aes(x = level, y = estimate)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), size = 0.3, width = 0.2) +
  scale_x_discrete(labels = c("None", "One", "Two", "Three")) +
  xlab('') + ylab('Marginal mean') +
  coord_flip() +
  theme_classic(base_size = 10) +
  scale_colour_grey() +
  easy_center_title() +
  ggtitle("(b) Disapproval ratings-based") +
  theme(
    axis.text.x = element_text(colour = "black"),
    axis.text.y = element_text(colour = "black"),
    axis.title.x = element_text(margin = margin(t = 10)),
    axis.ticks.y = element_blank()
  )

# Print Figure 5
pdf('figure5.pdf', width = 8.26, height = 5.82)
ggarrange(fig5_m1_plot, NULL, fig5_m2_plot, nrow = 1, common.legend = T, legend = "bottom", ncol=3, widths = c(1, 0.09, 1))
dev.off()

# =============================================================================
# Figure 6: Additive effect of multiple shared target and elite cue giver identities, conditional on perpetrator (partisanship)
# =============================================================================

# -----------------------------------------------------------------------------
# Calculate marginal means interactions
# -----------------------------------------------------------------------------

# Model 1
fig6_m1 <- cregg::cj(hr_survey, outcome1 ~ agent + type + scope + targ_nonstate + targ_race_match + targ_relig_match + targ_citiz_match + frame + elite_match + perp_targ_race + perp_targ_relig + perp_targ_citiz + perp_elite + shared_ident, id = ~id, estimate = "mm", by = ~perp_match)

# Model 2
fig6_m2 <- cregg::cj(hr_survey, outcome2 ~ agent + type + scope + targ_nonstate + targ_race_match + targ_relig_match + targ_citiz_match + frame + elite_match + perp_targ_race + perp_targ_relig + perp_targ_citiz + perp_elite + shared_ident, id = ~id, estimate = "mm", by = ~perp_match)

# Subset shared identity variable
group_vars <- c("shared_ident")
fig6_m1 <- fig6_m1[fig6_m1$feature %in% group_vars,]
fig6_m2 <- fig6_m2[fig6_m2$feature %in% group_vars,]

# Exclude NA values
fig6_m1 <- fig6_m1[!is.na(fig6_m1$estimate),]
fig6_m2 <- fig6_m2[!is.na(fig6_m2$estimate),]

# -----------------------------------------------------------------------------
# Prepare Figures
# -----------------------------------------------------------------------------

# Combine models
fig6_all <- rbind(fig6_m1, fig6_m2)

# Set factor order
fig6_all_perp_in <- filter(fig6_all, !grepl("perp_out", BY))
fig6_all_perp_in <- fig6_all_perp_in %>%
  mutate(level = factor(level, levels = c("shared_3", "shared_2", "shared_1", "shared_0")) %>% fct_rev())
fig6_all_perp_out <- filter(fig6_all, !grepl("perp_in", BY))
fig6_all_perp_out <- fig6_all_perp_out %>%
  mutate(level = factor(level, levels = c("shared_3", "shared_2", "shared_1", "shared_0")) %>% fct_rev())

# Separate models
fig6_m1_perp_in <- filter(fig6_all_perp_in, outcome == "outcome1")
fig6_m2_perp_in <- filter(fig6_all_perp_in, outcome == "outcome2")
fig6_m1_perp_out <- filter(fig6_all_perp_out, outcome == "outcome1")
fig6_m2_perp_out <- filter(fig6_all_perp_out, outcome == "outcome2")

# -----------------------------------------------------------------------------
# Create figures
# -----------------------------------------------------------------------------

# Create panel (a) Perpetrator in-group - disapproval forced-choice
fig6_m1_perp_in_plot <- ggplot(fig6_m1_perp_in, aes(x = level, y = estimate)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), size = 0.3, width = 0.2) +
  scale_x_discrete(labels = c("None", "One", "Two", "Three")) +
  xlab('') + ylab('Marginal mean') +
  coord_flip(ylim = c(0.38, 0.62)) +
  geom_hline(yintercept = 0.5, size = 0.2, color = "black") +
  theme_classic(base_size = 10) +
  easy_center_title() +
  ggtitle("(a) Disapproval forced-choice")

# Create panel (b) Perpetrator in-group - disapproval ratings-based
fig6_m2_perp_in_plot <- ggplot(fig6_m2_perp_in, aes(x = level, y = estimate)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), size = 0.3, width = 0.2) +
  scale_x_discrete(labels = c("None", "One", "Two", "Three")) +
  xlab('') + ylab('Marginal mean') +
  coord_flip(ylim = c(0.63, 0.72)) +
  theme_classic(base_size = 10) +
  easy_center_title() +
  ggtitle("(b) Disapproval ratings-based")

# Create panel (c) Perpetrator out-group - disapproval forced-choice
fig6_m1_perp_out_plot <- ggplot(fig6_m1_perp_out, aes(x = level, y = estimate)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), size = 0.3, width = 0.2) +
  scale_x_discrete(labels = c("None", "One", "Two", "Three")) +
  xlab('') + ylab('Marginal mean') +
  coord_flip(ylim = c(0.38, 0.62)) +
  geom_hline(yintercept = 0.5, size = 0.2, color = "black") +
  theme_classic(base_size = 10) +
  easy_center_title() +
  ggtitle("(c) Disapproval forced-choice")

# Create panel (d) Perpetrator out-group - disapproval ratings-based
fig6_m2_perp_out_plot <- ggplot(fig6_m2_perp_out, aes(x = level, y = estimate)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), size = 0.3, width = 0.2) +
  scale_x_discrete(labels = c("None", "One", "Two", "Three")) +
  xlab('') + ylab('Marginal mean') +
  coord_flip(ylim = c(0.63, 0.72)) +
  theme_classic(base_size = 10) +
  easy_center_title() +
  ggtitle("(d) Disapproval ratings-based")

# Print Figure 6
fig6_perp_in_plot <- ggarrange(fig6_m1_perp_in_plot, NULL, fig6_m2_perp_in_plot,
                               nrow = 1, common.legend = F, legend = NULL, ncol = 3, widths = c(1, 0.09, 1))
fig6_perp_in_plot<-annotate_figure(fig6_perp_in_plot, top = text_grob("Perpetrator in-group", size = 14))
fig6_perp_out_plot <- ggarrange(fig6_m1_perp_out_plot, NULL, fig6_m2_perp_out_plot,
                                nrow = 1, common.legend = T, legend = "bottom", ncol = 3, widths = c(1, 0.09, 1))
fig6_perp_out_plot<-annotate_figure(fig6_perp_out_plot, top = text_grob("Perpetrator out-group", size = 14))
pdf('figure6.pdf', width = 8.26, height = 5.82)
ggarrange(fig6_perp_in_plot, fig6_perp_out_plot,
          nrow = 2, common.legend = T, legend = "bottom", ncol=1)
dev.off()